#'
#' pairwise species association 
#' @aliases pairwise_association pairwise_association.community test_pairwise_association test_pairwise_association_community
#' 
#' @usage
#' pairwise_association(com,r,edgeCorr="WM")
#' test_pariwise_association(com_obs,com_rec,r=seq(0,30,0.5),method="MAD")
#' test_pairwise_association_community(com,nsim=99,r,dimyx=c(10,10),sigma=5)
#' onetoall_association(com,r,edgeCorr="WM")
#' test_pariwise_association_onetoall(com_obs,com_rec,r=seq(0,30,0.5))
#' 
#' 
#' @param com a community object
#' @param com_rec a list of reconstructed communities.
#' @param r the distance where the pariwise species association was evaluated, for trees, typical rmax is 25.
#' @param method which method, MAD or Gof test, should be used in the global significant test
#' @param edgeCorr edge correction method, only iso is the default setting. No other method implemented yet
#' @param focus_sp the character of focus species. if this is given, it will only calculate pairwise association involving this species.
#' @param nsim number of point pattern reconstructions used in testing the pairwise species association
#' @param dimyx,sigma the parameters used in kernel estiamte the density of individuals across the plot
#' @param alpha significant level used to calcualte the significance of the difference between observation and null model
#' @param tsigr vector of interested distance range from zero to do the Loosmore goodness of fit test
#'
#' @details
#' pairwise_association is used to calculate associations of all pair of species in the given community.
#' onetoall_association is used to calcualte association of a focal species with other neighbors.
#' 
#' test_pairwise_association is used to test the significance of each pairwise association by the given null communities.
#' 
#' test_pairwise_association_community is used to directly test the significance of each pairwise association.
#' it means this function will call pattern reconstruction to generate the null communities before test the
#' significance of pairwise species association.
#'
#' @return
#' A paas object
#' The test_pairwise_association returns a 'pairasso_com' object, which is a list contains a vector of all
#' pairwise species assocation by pcf, a vector of pairwise confidence interval and a vector of pvalues of 
#' pairwise species association and a vector of effect size adjusted by the normalization.
#'
#' The 'pairasso_com' object can be used in plot, sig_percent
#'
#'@examples
#'library(spatstat)
#'data(BCI)
#'
#'spab=species_abundance(BCI)
#'testsp=names(spab)[spab>500 & spab<610]
#'
#'subbci=get_populations(BCI,testsp)
#'
#'pairwise_association(subbci,r=seq(0,10,1),bw=1)
#'
#'splist = c("sp1", "sp2", "sp3")
#'abundance = c(200, 200,500)
#'obs_com=rcommunity(splist=splist,abundance=abundance)
#'rec_com=replicate(199,rcommunity(splist=splist,abundance=abundance),simplify=FALSE)
#'re=test_pariwise_association(obs_com,rec_com,method="Gof",mccores=1)
#'re2=test_pariwise_association(obs_com,rec_com,method="MAD",mccores=1)
#'
#'test_pariwise_association_onetoall(obs_com,rec_com)
#'
#'# a plot method for pairasso_com object
#'plot(re)
#'
#'#it needs a long running time, please wait......
#'test_pairwise_association_community(subbci,nsim=2,r=seq(0,10,1),tsigr=c(5,10,15))
#'
#'
#'

pairwise_association=function(data,...){
  UseMethod("pairwise_association")
}

pairwise_association.population=function(pop1,pop2,r){
  edgeCorr="isotropic"
  #step1 calculate observed cross pcf for each population
  ppp1=pop_to_ppp(pop1)

  #step2 return the cross pcf
}

# 
# pairwise_association.community=function(com,r,focus_sp=NULL,bw){
#   require(spatstat)
#   edgeCorr="isotropic"
#   #step1 calculate observed pcf for each pair of population
#   S=total_richness(com)
#   splist=species_list(com)
#   nr=length(r)
#   
#   com_ppp=com_to_ppp(com)
#   if(is.null(focus_sp)){
#     result=array(dim=c(S,S,nr))
#     spi=S+1
#   }else{
#     result=array(dim=c(S,nr))
#     spi=which(splist==focus_sp)
#   }
#   
#   for(i in 1:S){
#     if(!is.null(focus_sp) & i!=spi){
#       next()
#     }
#     for(j in 1:S){
#       if(i!=j){
#         if(is.null(focus_sp)){
#           result[i,j,]=pcfcross(com_ppp,splist[i],splist[j],r=r,correction=edgeCorr,bw=bw)$iso
#         }else{
#           result[j,]=pcfcross(com_ppp,splist[i],splist[j],r=r,correction=edgeCorr,bw=bw)$iso
#         }
#         
#       }
#     }
#   }
#   attr(result,"r")=r
#   return(result)
# }


#multicore verion
single_pair=function(i,ijsp,com_ppp,splist,r,edgeCorr,bw=NULL){
  pcfcross(com_ppp,splist[ijsp[i,1]],splist[ijsp[i,2]],r=r,correction=edgeCorr,bw=bw)$iso
}

pairwise_association.community=function(com,r,focus_sp=NULL,bw=NULL,mc.cores=1){
  require(spatstat)
  edgeCorr="isotropic"
  #step1 calculate observed pcf for each pair of population
  S=total_richness(com)
  splist=species_list(com)
  nr=length(r)
  
  com_ppp=com_to_ppp(com)
  if(is.null(focus_sp)){
    result=array(dim=c(S,S,nr))
    spi=S+1
    ijsp=expand.grid(i=1:S,j=1:S)
    
  }else{
    result=array(dim=c(S,nr))
    spi=which(splist==focus_sp)
    ijsp=expand.grid(i=spi,j=1:S)
  }
  
  ijsp=ijsp[-which(ijsp$i==ijsp$j),]
  
  if(mc.cores==1){
    all_re=lapply(1:(dim(ijsp)[1]),single_pair,ijsp=ijsp,com_ppp=com_ppp,
                  splist=splist,r=r,edgeCorr=edgeCorr,bw=bw)
  }else{
    all_re=mclapply(1:(dim(ijsp)[1]),single_pair,ijsp=ijsp,com_ppp=com_ppp,
                  splist=splist,r=r,edgeCorr=edgeCorr,bw=bw,mc.cores = mc.cores)
  }
  
 
  for(i in 1:(dim(ijsp)[1]) ){
      if(is.null(focus_sp)){
        result[ijsp[i,1],ijsp[i,2],]=all_re[[i]]
      }else{
        result[ijsp[i,2],]=all_re[[i]]
      }
  }
  
  attr(result,"r")=r
  return(result)
}



onetoall_association=function(com,r,focus_sp=NULL){
  require(spatstat)
  edgeCorr="isotropic"
  #step1 calculate observed pcf for each pair of population
  S=total_richness(com)
  splist=species_list(com)
  nr=length(r)
  
  com_ppp=com_to_ppp(com)
  if(is.null(focus_sp)){
    result=array(dim=c(S,nr))
  }else{
    result=array(dim=c(1,nr))
  }
  
  
  if(is.null(focus_sp)){
    for(i in 1:S){
      isfocal=com_ppp$marks$species==splist[i]
      I=which(isfocal)
      J=which(!isfocal)
      result[i,]=pcfmulti(com_ppp,I,J,r=r,correction=edgeCorr)$iso
    }
  }else{
    isfocal=com_ppp$marks$species==focus_sp
    I=which(isfocal)
    J=which(!isfocal)
    result[1,]=pcfmulti(com_ppp,I,J,r=r,correction=edgeCorr)$iso
  }
  
  
  attr(result,"r")=r
  return(result)
}

#test the significance of the association between a focal species and other neibors
test_pariwise_association_onetoall=function(com_obs,com_rec,r=seq(0,30,0.5),method="MAD"){
  splist=species_list(com_obs)
  S=total_richness(com_obs)
  nsim=length(com_rec)
  nr=length(r)
  #calcualte the observed pairwise assocaition 
  obs_result=onetoall_association(com_obs,r=r)
  
  #save for effect size
  ef_result=array(NA,dim=dim(obs_result))
  #stand devation 
  sd_null=ef_result
  
  #calcualte the pairwise association with random community
  null_result=array(dim=c(nsim,nr))
  result_pvalue=matrix(nrow=S,ncol=1)
  result_conf=array(dim=c(S,2,nr))
  
  for(i in 1:S){
    focus_sp=splist[i]
    sel_obs=com_obs$species==focus_sp
    for(j in 1:nsim){
      one_rec_com=com_obs
      sel_sim=com_rec[[j]]$species==focus_sp
      one_rec_com$x[sel_obs]=com_rec[[j]]$x[sel_sim]
      one_rec_com$y[sel_obs]=com_rec[[j]]$y[sel_sim]
      null_result[j,]=onetoall_association(one_rec_com,r=r,focus_sp=focus_sp)
    }
    
    #calcualte the confidence interval
    result_conf[i,,]=apply(null_result,2,conf)
    sd_null[i,]=apply(null_result,2,sd,na.rm=TRUE)
    #calculate the effect size based on simulations under the null communities
    ef_result[i,]=(obs_result[i,]-apply(null_result,2,mean,na.rm=TRUE))/sd_null[i,]
    
    #do the MAD test or Loosmore goodness of fit test
    if(method=="Gof"){
      result_pvalue[i,]=lmGofTest(cbind(obs_result[i,],t(null_result)))
    }else{
      re_sigtest=MAD(cbind(obs_result[i,],t(null_result)),alpha=0.05,retunfpvalue=TRUE,obsadj=TRUE)
      ef_result[i,]=attr(re_sigtest,"obsadj")
      result_pvalue[i,]=as.numeric(re_sigtest)
    }
    
    
  }
  
  all_result=list(obs=obs_result,conf=result_conf,pvalue=result_pvalue,effectsize=ef_result,sdnull=sd_null)
  class(all_result)=c("onetoallasso_com")
  return(all_result)
  
}



#test the significance of the pairwise association by many given reconstructed community
test_pariwise_association=function(com_obs,com_rec,r=seq(0,30,0.5),bw=0.5,mccores=40,
                                   method="MAD",obs_result=NULL,mc=FALSE,two_side=TRUE){
  splist=species_list(com_obs)
  S=total_richness(com_obs)
  nsim=length(com_rec)
  nr=length(r)
  #calcualte the observed pairwise assocaition 
  if(is.null(obs_result)){
    obs_result=pairwise_association(com_obs,r=r,bw=bw,mc.cores=mccores)
  }
  
  #save for effect size
  ef_result=array(NA,dim=dim(obs_result))
  sd_null=ef_result
  norm_p=ef_result
  
  #calcualte the pairwise association with random community
  null_result=array(dim=c(nsim,S,nr))
  result_pvalue=matrix(nrow=S,ncol=S)
  result_conf=array(dim=c(S,S,2,nr))
  
  for(i in 1:S){
    focus_sp=splist[i]
    sel_obs=com_obs$species==focus_sp
    com_sel=com_obs[sel_obs,]
    #browser()
    if(.Platform$OS.type=="windows" | !mc){
      tempre=lapply(com_rec,foronesp_inner,focus_sp=focus_sp,com_sel=com_sel,r=r,bw=bw)
    }else{
      require(parallel)
      #adjusted for multicore calculation, it will be useful for many reconstructions
      tempre=mclapply(com_rec,foronesp_inner,focus_sp=focus_sp,com_sel=com_sel,r=r,bw=bw,
                                     mc.cores=mccores)
    }
    for(j in 1:nsim){
      null_result[j,,]=tempre[[j]]
    }
    
    for(k in 1:S){
      if(k!=i){
        #calcualte the confidence interval
        result_conf[i,k,,]=apply(null_result[,k,],2,conf)
        sd_null[i,k,]=apply(null_result[,k,],2,sd,na.rm=TRUE)
        #calculate the effect size based on simulations under the null communities
        ef_result[i,k,]=(obs_result[i,k,]-apply(null_result[,k,],2,mean,na.rm=TRUE))/sd_null[i,k,]
        #print(c(i,k))
        #if(i==1 & k==15) browser()
        #normality test
        norm_p[i,k,]=c(NA,apply(null_result[,k,-1],2,function(x) ifelse(all(x==x[3]),NA,shapiro.test(x)$p.value)))
        
        #do the Loosmore goodness of fit test
        #result_pvalue[i,k]=lmGofTest(cbind(obs_result[i,k,],t(null_result[,k,])))
        if(method=="Gof"){
          result_pvalue[i,k]=lmGofTest(cbind(obs_result[i,k,],t(null_result[,k,])))
        }else{
          re_sigtest=MAD(cbind(obs_result[i,k,],t(null_result[,k,])),alpha=0.05,two_side=two_side,
                         retunfpvalue=TRUE,obsadj=TRUE)
          ef_result[i,k,]=attr(re_sigtest,"obsadj")
          result_pvalue[i,k]=as.numeric(re_sigtest)
        }
      }
    }
    
  }
  
  all_result=list(obs=obs_result,conf=result_conf,pvalue=result_pvalue,effectsize=ef_result,sdnull=sd_null,norm_p=norm_p)
  class(all_result)=c("pairasso_com")
  return(all_result)
  
}

#extact from the test_parise_association for multicore calculation
foronesp_inner <- function (one_rec_com,focus_sp,com_sel,r,bw) {
  #one_rec_com=com_rec[[j]]
  #browser()
  #one_rec_com[one_rec_com$species==focus_sp,]=com_sel
  one_rec_com=rbind(one_rec_com[one_rec_com$species!=focus_sp,],com_sel)
  tp=pairwise_association(one_rec_com,r=r,focus_sp=focus_sp,bw=bw)
  return(tp)
}


plot.pairasso_com=function(re,type="pair_conf",pair=NULL){
  S=dim(re$pvalue)[1]
  r=attr(re$obs,"r")
  if(type=="pair_conf"){
    #draw all pair association plots
    if(is.null(pair)){
      par(mfrow=c(S,S))
      for(i in 1:S){
        for(j in 1:S){
          if(i==j){
            plot(x=1,y=1,type="n")
            next()
          }
          confj=re$conf[i,j,,]
          plot(x=range(r),y=range(confj[,-1],na.rm=T),type="n",xlab="r",ylab="pcf")
          lines(x=r,y=re$obs[i,j,])
          lines(x=r,y=confj[1,],col="grey")
          lines(x=r,y=confj[2,],col="grey")
        }
      }
      
    }
  }else if(type=="perc_sig"){
    perc_up=perc_down=rep(0,length(r))
    
    for(i in 1:S){
      for(j in 1:S){
        if(i!=j){
          confj=re$conf[i,j,,]
          perc_up=perc_up+re$obs[i,j,]>confj[2,]
          perc_down=perc_down+re$obs[i,j,]<confj[1,]
        }
      }
    }
    perc_up=perc_up/S
    perc_down=perc_down/S
    plot(x=range(r),y=c(0,1),type="n",xlab="r",ylab="Percentage")
    lines(x=r,y=perc_up,type="o",pch=1)
    lines(x=r,y=perc_down,type="o",pch=2)
    lines(x=r,y=1-perc_down-perc_up,type="o",pch=3)
    return(data.frame(r=r,sig_up=perc_up,insig=1-perc_down-perc_up,sig_down=perc_down))
  }
}

sig_percent=function(re,siglevel=0.05,adjust=FALSE){
  n=dim(re$pvalue)[1]
  return(sum(re$pvalue<=siglevel,na.rm=TRUE)/(n*(n-1)))
}

#test the significance of the association by point pattern reconstruction
test_pairwise_association_community=function(com,nsim=99,r=seq(0,30,0.5),dimyx=c(10,10),sigma=5,alpha=0.05,tsigr=seq(5,30,5)){
  S=total_richness(com)
  splist=species_list(com)
  result=list()
  cat(paste("Total",S,"species,","finished species:"))
  for(i in 1:S){
    result[[i]]=test_pairwise_association_onesp(splist[i],com,nsim,r,dimyx,sigma,alpha,tsigr)
    cat(paste(i,",",sep=""))
  }
  names(result)=splist
  return(result)
}

#test the significance of the association by point pattern reconstruction
test_pairwise_association_onesp=function(spname,com,nsim=99,r,dimyx=c(10,10),sigma=5,alpha=0.05,tsigr=seq(5,30,5)){
  pop=get_population(com,spname)
  #reconstruct nr number of a population
  rec_pops=ppreconstruction(pop,nrc=nsim,intmap=TRUE,keepfiles=FALSE,dimyx=dimyx,sigma=sigma)
  #calcualte the pariwse association of this observed and reconstructed populations with other populations
  S=total_richness(com)
  splist=species_list(com)
  nr=length(r)
  
  com_ppp=com_to_ppp(com)
  sel=com_ppp$marks$species==spname
  spi=which(spname==splist)
  result=array(dim=c(S,nr,nsim+1))
  
  #simulated value
  for(i in 1:(nsim+1)){
    if(i==1){
      com_ppp$x[sel]=pop$x
      com_ppp$y[sel]=pop$y
    }else{
      com_ppp$x[sel]=rec_pops[[i-1]]$x
      com_ppp$y[sel]=rec_pops[[i-1]]$y
    }
    for(j in 1:S){
      if(j!=spi){
        #the first i is the observed value
        result[j,,i]=pcfcross(com_ppp,spname,splist[j],r=r,correction="isotropic")$iso
      }
    }
  }
  
  re=list()
 
  pvalues=matrix(nrow=S,ncol=length(tsigr))
  colnames(pvalues)=tsigr
  rownames(pvalues)=splist
  
  if(max(r)<max(tsigr)){
    tsigr=tsigr[tsigr<=max(r)]
    warnings("max value of tsigr larger than max r")
  }
  
  for(i in 1:S){
    if(i!=spi){
      #calcualt the confidence interval   
      confi=apply(result[i,,-1],1,conf,alpha=alpha,na.rm=TRUE)
      re[[i]]=rbind(confi,result[i,,1]) #obs, min, max
      rownames(re[[i]])=c("obs","min","max")
      colnames(re[[i]])=r
      #Loosmore goodness of fit test for all tsigr
      for(j in 1:length(tsigr)){
        selr= r<= tsigr[j]
        pvalues[i,j]=lmGofTest(result[i,selr,])
      } 
    }else{
      re[[i]]=NA
    }
  }
  
  names(re)=splist
  #save these pairwise associations
  attr(re,"pvalues")=pvalues
  return(re)
}


#plot the pairwise association
plot.paas=function(pairass){
  
}
